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ABSTRACT 

Context. The development of parallel supercomputers allows today the detailed study of the collapse and the fragmentation of prestel- 
lar cores with increasingly accurate numerical simulations. Thanks to the advances in sub-millimeter observations, a wide range of 
observed initial conditions enable us to study the different modes of low-mass star formation. The challenge for the simulations is to 
reproduce the observational results. 

Aims. Two main numerical methods, namely AMR and SPH, are widely used to simulate the collapse and the fragmentation of 
prestellar cores. We compare thoroughly these two methods within their standard framework. 

Methods. We use the AMR code RAMSES and the SPH code DRAGON. Our physical model is as simple as possible, and consists of 
an isothermal sphere rotating around the z-axis. We first study the conservation of angular momentum as a function of the resolution. 
Then, we explore a wide range of simulation parameters to study the fragmentation of prestellar cores. 

Results. There seems to be a convergence between the two methods, provided resolution in each case is sufficient. Resolution criteria 
adapted to our physical cases, in terms of resolution per Jeans mass, for an accurate description of the formation of protostellar cores 
are deduced from the present study. This convergence is encouraging for future work in simulations of low-mass star formation, 
providing the aforementioned criteria are fulfilled. 

Key words. Stars: formation - Methods : numerical - hydrodynamics 



1. Introduction Nowadays, two completely different numerical method are used 

with sufficient accuracy: 

. . Star formation is known for being the place of extreme variations 

_ > in length and density scales. Although it is established that stars ^ ^^j^. Adaptive Mesh Refinement method for Eulerian grids 

form in dense cores, the non-linear evolution makes it difficult 2. SPH: Smoothed Particle Hydrodynamics method for a 

to perform accurate calculations of the collapse and the frag- Lagrangian approach 
mentation of a prestellar core. The star formation process is the 

outcome of complex gas dynamics involving non-linear interac- systematic comparison between the two methods has 

tions of gravity, turbulence rnagnetic feW an^ra^^ ^^^^ ^^^^ ^.^^ low-mass star formation calculations. However, 
theo r etical pioneer works by ILarsonl (119691) . IPenstonI 019691) or 



77, — I nrr^* . ' ^r; ~ — T"; — , ~, — a lot of numerical works have been carried out and some of 

\Shd are among the many illustration s of the hig hcgm- ^^^^ ^^^^^^ ^^^^ calculations for convergence testing and 

plexity of the gravitational collapse. Recently, Klein et al. ( 2007') i^te^code c omparisons. The mos t famous model was first cal- 

pointed out that developing a theory tor low-mass star formation 1 » ju Fd~^ — p~dj — u~- — n hnid. a ■ »u ■» u 

^ . ^ , ?. ,. , r , culated bv B oss & Bodenheimerl 019791) and since then, it has 

remains one of the most elusive and important goals of theoret- ^^^^ recalculated by several auth ors w ith even higher spa- 

ical astrophysics. The computational challenge stems from the tial resolution (e.g 'Bate & Burkert" 1997; Truelove et al.lj[l99i 

fact that star formation occurs in clouds over many orders of Ixitsionas & Whitworth 2002; Arreaga-Gaixia et al.. .20071) ^ 
magiiitude in spatial and density scales. Following the gravita- j^^^ generated a lot of detailed investigations 

tional collapse while resolving precisely the Jeans length, which ^^^3^^^^^^ ^^^^^^ .^^^^ neighbors 

scales as A, oc p-i 2 for an isothermal gas, is a major difficulty (^ombardi et al. 1999; Rasio 1999; Attwood et al. 2007), and 

tor numerical simulations. criteria for n umerical convergence have been extracted from 

Different approaches are used to study star formation t^ese studies. iNelsonI (12006 ) performed a large investigation 
through numerical simulations and include more and more de- influence of fliese parameters on disk fragmentation, 

tailed physics. One key question resides in the validation of ^nd concluded that the better the resolution the later the frag- 

the numerical methods used to study low-mass star formation, mentation. Dehnen (2001) investigated the optimal gravitational 

force softening necessary in three-dimensional A^-body codes. 

Send offprint requests to: B. Commer§on iBate & BurkertI ([T997i) provide a minimum resolution criterion 
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for SPH calculations with self-gravity to accurately model frag- 
mentation. Less studies have been performed with AM R since 
AMR codes have become available recently. Truelov e et alj 
(Il997h give an empirical criterion for the Jeans length resolution 
in AMR calculations to avoid spurious numerical fragmentation. 

There are not much of direct comparison between SPH and 
AMR calculations. Comparison in the context of cosmological 
simulations has been done throiigh the Santa Barbara Cluster 
Comparison Project (iFrenk et al.l Il999l) . iFromang et al.l (l2006l) 
compares quite successfully AMR hydrodynamical c ollaps e 
calculations with the ones of iHosking & WhitworthI (|2004|) . 
using the SPH method. 

In the present paper, we compare thoroughly the two ap- 
proaches in the context of low-mass prestellar core formation. 
We stress that the main goal of this paper is to investigate 
whether convergence can be achieved between the two methods. 
We have conducted calculations over a wide range of numerical 
resolution parameters, in order to study the dependency of angu- 
lar momentum conservation and fragmentation on physical and 
numerical initial conditions. We then derive resolution criteria 
necessary to describe accurately prestellar core formation. 

The paper is organized as follows: in §2 we briefly introduce 
our collapse model. In §3, we present the two codes used for our 
comparative study as well as our initial numerical conditions 
and the criteria we fulfill to resolve gravitational collapse. The 
problem of angular momentum conservation is examined in 
detail §4. In §5, we tackle the fragmentation issue and explore 
the dependency of the results on the numerical parameters. First, 
we study the numerical convergence of AMR and SPH calcu- 
lations separately. Then, we compare the respective converged 
calculations. This convergence study is done for different test 
cases. In section 6 we conclude this paper by deriving for each 
method empirical required numerical criteria for an accurate 
description of gravitational collapse and fragmentation adapted 
to our test cases. 

The convention in this paper is to call "particles" the SPH 
particles and "cells" the AMR cells in order to avoid confusion. 

2. Definitions of the test cases 

2.1. Model 

To make comparison between codes easier, we adopt simple ini- 
tial condit ions, simil ar to those chosen in previo us studies (e.g. 
iBoss & Bodenheimerll979HBate & Burker3[T997l) . We consider 
an uniform-density sphere of molecular gas of initial radius Rq, 
rotating around the z-axis with an uniform angular velocity Qq, 
in order to minimize the loss of angular momentum due to fric- 
tion. We fix the cloud mass at Mq - 1 Mq and the tempera- 
ture at 10 K. For a mixture of molecular hydrogen, helium and 
heavy elements, this corresponds to an isothermal sound speed 
of Co ~ 0.19 km.s"'. For the case where fragmentation occurs, 
we use a m = 2 azimuthal density perturbation. 

The initial energy balance of our model is determined by two 
dimensionless parameters corresponding to the ratio between the 
thermal energy and the gravitational energy 

5RokT 
2 GMoju;«h 

and to the ratio of the rotational and the gravitational energy 
1 Rpl 



Since we use a constant initial mass of 1 Mq and a con- 
stant temperature, changing one of the two parameters, namely 
a, gives the sphere radius Rq. The higher a, the larger Rq. The 
angular velocity is given by the parameter jS. 

2.2. The barotropic equation of state 

In order to mimic the thermal behaviou r of a star-forming gas, 
we use a ba rotropic equation of sta t e (cf.lBonnellll994l) . lTohlind 
(1198 2') and 'Masu naga & Inutsukal ilOOU) showed that the core 
follows closely a barotropic equation of state, providing a good 
approximation without resolving radiative transfer. We use 



p 



2/3 



(3) 



where is the sound speed and pc = 10 g.cm" is the critical 
density which corresponds to the transition from an isothermal 
to an adiabatic state (Larson 1969). 

At low densities, p p^, Cs ~ Co = 0.19 km.s"'. The 
molecular gas is able to radiate freely by coupling thermally 
to the dust and therefore remains isothermal at 10 K. At high 
densities p > p^, we assume that the cooling due to radiative 
transfer is trapped by the dust opacity. Therefore, P cc p^l^ 
which corresponds to an adiabatic monoatomic gas with adi- 
abatic exponent y - 5/3. Note that molecular hydrogen be- 
haves like a monoatomic gas until the temperature reaches sev- 
eral hundred Kelvin, since the rotational degrees of freedom 
are not excited at lower temperatu res, and hence y = 5/3 is 
the appropriate adiabatic ex ponent dWhitworth & Clark3ll997l; 
iMasunaga & Inutsukall2000h . 

3. Numerical methods and initial conditions 

3.1. AMR: Adaptive Mesli Refinement 
3.1.1. A brief history 

The Adaptive Mesh Refinement method is one of the most 
promising numerical methods to so lve the fluid equations . 
The technique was first introduced in iBerger & Oligeii (|1984|) . 
Originally, the AMR method was an Eulerian hydrodynamical 
scheme, with a hierarchy of nested grids covering high res- 
olution regions of the flow. This first AMR structure, called 
"patch-based AMR", consists of building blocks of the compu- 
tational grid as rectangular patche s of various size s. An alter- 
native method was proposed (i.e. lKhokolovlll998l) . the "tree- 
based" AMR, where the parent cells are refined into chil- 
dren cells on a cell-by-cell basis. These adaptive mesh struc- 
tures are coupled with grid-based fluid dynamics schemes han- 
dling high-resolution shock capturing. Nowadays, high order 
Godunov methods appear to be amongst the best sc hemes to 



capture discontinuities within only a few cells (e.g. T evssiei 
2002; Matsumoto & Hanawa2003:.Ziegler.2005: Fromang et al 



12006^ . 



3.1.2. Tine RAMSES code 



3 GMo 



(2) 



In th is paper, we use the AMR code RAMSES dTevssied 
I2002h . which integrates the "tree-based" data structure allowing 
recursive grid refinements. RAMSES uses a second order 
Godunov hydrodynamical scheme coupled with a gravity solver 
Furthermore, it has the possibility to use variable timesteps at 
each refinement level. Concerning time integration, RAMSES 



B. Commerjon et al.: Protostellar collapse: A comparison between SPH and AMR calculations. 



3 



uses a second-order midpoint scheme, where positions and 
velocities are updated by a predictor-corrector step. Recently, 
an ideal MHD versi on of RAMSES has been developed by 
lFromangetaLl(l2006l) . 

The Godunov hydrodynamical solver is able to capture dis- 
continuities with a high precision level. The equations solved in 
RAMSES are the Euler equations in their conservative form 



dt 



[pu] = 0, 



, u -I- PI] = -pV<D, 



— .Vbu( 
dE 

— -H V [u (E -H P)] = -pu ■ V(D. 

dt 



(4) 



(5) 



(6) 



where P is the gas pressure (I is the identity matrix), p the den- 
sity, u the velocity, E the total energy density and O the grav- 
itational potential. The system of equations is closed with the 
barotropic equation of state ([3]). 

One of the main advantages of solving the Euler equations in 
their conservative form is that no energy sink due to numerical 
errors can alter the flow dynamics (ignoring the source terms due 
to gravity). Equations |4l|5] and |6] are solved with a Lax-Friedrich 
Riemann solver, known to be one of the most simple and robust 
scheme. In Appendix ICl we report on the influence of using the 
Roe solver which is less diffusive but considerably more com- 
plex and numerically expensive. 

The timestep is determined independently for each refine- 
ment level Hi, using standard stability constraints for the hy- 
drodynamical solver Each level £i evolves according to its own 
timestep. 

3.1 .3. Initial conditions for RAMSES 

One practical limitation of AMR codes is the use of Cartesian 
grids. RAMSES works with a cubic volume, so that a part of 
the calculation box is lost when we describe a sphere. The outer 
region of the sphere is also at a 10 K temperature but is 100 
times less dense. Therefore, the outer gas has no effect on the 
dynamics of the sphere since the two parts are well separated. 
The sphere radius is equal to a quarter of the box length in order 
to minimize border effects. 



3.2. SPH: Smoothed Particles Hydrodynamics 
3.2.1 . A brief overview 

SPH is the most popular fully Lagrangian method used to de- 
scribe gravitational collapse because of its simplicity for 3D 
codes and its versatility to incorporate self-gravity. SPH was first 
designed to sim ulate nonaxisymmetric phenomena for astro- 
physical gases (,Lucvlll977t iGingold & MonaghanI [19771) . This 
method is easy to work with and can give rapidly reasonably 
accurate results. SPH is economic in handling hydrodynamical 
flows that have near empty regions. It does not need a grid to 
calculate spatial derivatives, but consists of a set of discrete par- 
ticles describing the state of the fluid. The spatial derivatives 
are found by analytical differentiation of interpolation formulae. 
SPH particle / should not be perceived as a real fluid element, but 
as a mathematical entity with coordinates r,-, velocity v,-, mass m,- 
(i.e. m since all particles have the same mass in the present cal- 
culations) and thermal energy e,. The evolution of the fluid is 



determined by following the motion of the particles, under the 
influence of interparticle forces which represent the effects of 
pressure, viscosity (see below) and self-gravity. 

The main advantage of SPH is its strict Galilean-invariant 
property and its simplicity. Resolution elements are then con- 
centrated in high density regions in SPH methods . The stan- 
dard SPH formalism uses artificial viscosity for the hydrody- 
namics. Some altern ative formalism such as Godunov SPH has 
been proposed (e.g. Ilnutsukalll994l) in order to avoid the use 
of artificial viscosity, but these methods are not yet mature. 
SPH has been used by several authors to study fragmenta- 
tion(Bonnell 1994; Bat e & Burkert|[T997HGoodwin et al.ll2004l: 
Hennebelle et aL,2004i) . 

3.2.2. Hydrodynamical method for DRAGON 

We use the standard SPH code DRAGON (iTurner et al.lll995l: 
iGoodwin et alj|200 4). i.e in its most simple version. In standard 
SPH, the integral interpolant for the variable A(r,) is approxi- 
mated by a summation interpolation over the particle's nearest 
neighbors: 



Air,) 
V V(r>) 



(7) 



where A(ry) is the value associated with particle j, hij - {hi + 
hj)/2 and hi is the adaptive smoothing length of particle /, de- 
fined such that the particle kernel volume contains a constant 
mass, i.e. a constant number of neighbors A^n- The interpolation 
mass is then given by mNt^. 

We use a standard artificial viscosity scheme 
(IGingold & Monaghanlll983l) : 



with 



^ f(-fl 1 CojjMij + aififj) Ipu if Mij 



ru < 0, 



otherwise. 



hiju^ 



(8) 



(9) 



where u,; = U;-Uys r,y - r,-ry and Cqjj and pij denote arithmetic 
means of the isothermal sound speed and density of the particles 
/ and / The free parameters ai, and regulate the strength 
of the viscosity. In our case, we have the combination oi = 1, 
a2 - 2 and = 0.1/;,/. This artificial viscosity is the subject of 
several discussions where authors suggest alternative artificial 
viscosity. One common improvement is to use a time-dependent 
viscosity (iMorris & Monaghanlll997h . We look at the influence 
of the viscosity scheme in AppendixiBl 

The gravitational force between a pair of particles obeys 
a simple inverse-square law, unless the particles are very 
close. Under this circumstance, the gravity force has to be 
softened to avoid vio lent two-body interaction. According to 
iBate & Burkertl (Il997h . the gravitational softening length should 
be the same as the hydrodynamic smoothing length. Calculation 
of the gravitational acceleration of an SPH particle is speeded 
up using an octal Spatial Tesselation Tree (STT) (iHernquisj 
1987) and accounts for the quadrupole moments of the mass 
distributions. DRAGON benefits from the implementation of 
"sink particles" creation (iBate et al.l 11995 ^. used to continue 
the calculations without resolving processes on extremely short 
time-scales. Last but not least, DRAGON uses multiple-particle 
timesteps. 



4 



B. Commerjon et al.: Protostellar collapse: A comparison between SPH and AMR calculations. 



Finally, we allow a variation of the number of neighbors AA^n 
less than 10% of A^n in our SPH calculations, i.e AAn - 5 when 
An - 50. lAttwood et all (|2007|) shows that the smaller AAn, the 
less diffusive is SPH. Ideally, AAn should be set to 0. We report 
on the influence of setting AAn = for local angular momentum 
conservation in AppendixiBl 

3.2.3. Initial conditions for DRAGON 



We us e the method originally presented by IWhitworth et al.l 
(Il995h to obtain an initial particle distribution which consist in 
settling the positions of the particles randomly settled by using 
only hydrodynamical forces over a few timesteps. Another ap- 
proach sometimes used is to take initial hexagonal-close-paced 
lattice of SPH particles to generate initial conditions. However, 
standard SPH calculations st art from noisy initial particle dis- 
tributions at present time (e.g. lArreaga-Garcia et al.i2007 :). Note 
that we do not need intercloud and external particles to confine 
the ones within the sphere since our model is initially far from 
equiUbrium. 

3.3. The Jeans criterion in numerical codes 

3.3.1 . Refinement criterion for the AMR method 

Our refinement criterion is based on the Jeans length resolution 
which is necessary to treat accurately gravitational collapse. We 
impose a minimum number of points Aj per Jeans length Aj. The 
cells' dimensions must be smaller than a constant fraction of the 
local Jeans length. The dimension of cells belonging to the ^, 
refinement level is Lbox/2'^', where Lbox is the physical length of 
the simulation box. The mesh is locally refined in order to satisfy 
the local Jeans criterion: 

^box ^] 



(10) 



iTruelove et al.l d 199 71) defined a minimum resolution condi- 
tion for the validity of grid-based simulations aimed at modeling 
the collapse of a molecular cloud core, namely Aj > 4. This 
condition ensures that the collapse is of physical rather than of 
numerical origin. 

3.3.2. Jeans length description with a SPH code 

In standard SPH, the resolution in mass is fixed and thus the 
Jeans length resolution deteriorates with increasing density for 
an isothermal gas. The minimum reso lvable mass mu s t then 
be larger than the interpolation mass. iBate & Burkeri (Il997h 
showed that the behaviour of a Jeans-mass clump of gas with 
radius ~ h is dominated by the numerical implementation. 
iBate & Burk ert (1997) take the smallest mass that can be re- 
solved in SPH calculations to be equal to the mass of ~ 2An 
particles. According to this criterion, we can determine an initial 
number of SPH particles necessary to solve the Jeans length in 
the simulations. 

The Jeans mass is Mj ~ 6G^^^^p^^^^Cl and the minimum 
resolvable mass is Mres - wAn. Hence, we can define a Jeans 
condition corresponding to the minimum value of C^p^^^^, given 
by the barotropic equation of state (O, i.e. l^^^C^Pc^^^: 



m < m„ 



5.35 X 10- 



■M(7 



0- 



(11) 



2ANG3/2py' 

Considering an initial spherical mass Mq - 1 M0, the initial 
number of particles Ap has to satisfy Ap > Mo/m^n-^j^ ~ 9300 if 



An = 50. This is the critical number of particles used in SPH cal- 
culations to study the collapse of a dense core. We have in that 
case exactly 2A n (i.e. two resolution elements) particles per crit- 
ical Jeans mass. lHubber et al1(l2006l) shows that with this numer- 
ical resolution, standard SPH will capture fragmentation which 
is genuine and resolved. 

As mentioned before, the mass resolution is fixed in stan- 
dard SPH. In term of Jeans mass, the resolution is therefore high 
at the beginning of the simulation and decreases when the den- 
sity increases up to the critical density. It is nevertheless instruc- 
tive to have a means of comparing the SPH and AMR resolu- 
tion. We therefore define for the SPH the parameter Aj such 
that Aj - M] I Mres is the number of resolution element per 
Jeans mass. This number is computed at the critical density p^, 
i.e. the most unfavorable case for the SPH. In the dense core 
where p > pc, the parameter Aj enables us to compare the res- 
olution achieved by SPH and A MR, respectively. Us ing particle 
splitting refinement in SPH (e.g. iKitsionas & WTiitworth .20021) 
would improve its resolution. Indeed, particle splitting is an eco- 
nomic way to increase the local resolution and thus to avoid vi- 
olating the Jeans condition in collapse simulations. However, as 
mentioned in the introduction, the aim of the present paper is to 
compare the AMR and the SPH within their standard implemen- 
tation. 



4. Free-fall time and angular momentum 
conservation 

We start by comparing the global properties of the collapse in 
the two codes in the simple case of an uniform-density sphere 
collapse with no perturbation. We look at the collapse time, the 
accretion shock and finally the angular momentum conservation. 

We carried out a first set of simulations within a wide range 
of resolution parameters. The initial sphere is set up by parame- 
ters a = 0.65 corresponding to an initial radius Rq = 9.2 x 10'^ 
cm and a density po ~ 6.02 x 10"'^ g.cm""'. The corresponding 
free-faU time is tff = (3;r/32Gpo)'^^ ~ 86 kyr. 



Table 1. Summary of the different simulations for the case with 
no rotation (left table: SPH; right table: AMR). Ai for the AMR 
calculations gives us the number of cells describing the initial 
sphere. 



5 X 10^' 50" 

1 X lO"* 50 
5 X 10^ 50 

2 X 10^ 50 
5 X 10' 50 



A'j to (kyr) 



1.86 
2.33 
4. 

6.35 
8.61 



108.6 
102 
94.1 
91.5 
90.6 



f . 


N, 


A'j 


to (kyr) 


5 


2 145 


10 


109 


6 


17 160 


10 


98 


7 


137 260 


10 


95 



4. 1 . Free-fall time 

The first step is to compare the calculations collapse time when 
the initial sphere is not rotating (i.e. /3 = 0). Note that since 
a is large and since we use a barotropic equation of state, we 
expect to find a value larger than the free-fall time. Then, we 
use as a reference time to ^ ( tff) for which p^ax = Pc- Table 
[1] gives collapse times to as a function of resolution parameters 
for AMR and SPH calculations. AMR calculations have been 
run with Aj = 10 and ^min = 5,6 and 7, SPH calculations with 
Ap ranging from 5 x 10^ to 2 x 10^ and a number of neighbors 
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log(r 



Fig. 1. (a): Density profiles at to as a function of the radius for the 
case with /3 - 0.01 and for the two most resolved simulations, 
namely (^ain - N] - 10 for the AMR (squares) and A^p - 
5 X 10^ A^N = 50 for the SPH (crosses). FigUlb): illustration 
of the accretion shock: radial velocity profiles for the two most 
resolved AMR and SPH simulations in the xy-plane at ~to + 1 .2 
kyr. 



A^N - 50. We note that with Ap = 5 X 10^, we do not satisfy the 
lBate&Burker3 (119971) criterion, but the mass of our resolution 
element, i.e. the sphere containing the An neighbors, is smaller 
that the critical Jeans mass. 

With increasing numerical resolution, the numerical time to 
decreases and seems to converge toward a value slightly greater 
than the free-fall time. Time to changes by less than 5% between 
AMR calculations with fn,in - 6 and 7 and SPH calculations 
with Ap = 5 X 10'* and 2 X 10^. Dynamical times to in SPH 
calculations are closer to the free-fall time than the AMR ones. 
This is partly due to the higher initial resolution in SPH. 

4.1 .1 . Collapse and accretion shock with rotation /3 = 0.01 

The gas sphere is now in solid rotation around the z-axis. We set 
y6 = 0.01, corresponding to an orbital time trot = 2.8 x lO-' kyr. 
Table|2]summarizes the different SPH and AMR calculations run 
for this case. In order to illustrate the core resolution, we give 
the quantity Acore representing the number of cells/particles with 
density p > 1 x 10"''' g.cm"-' at to. AMR simulations have been 
performed with different minimum refinement levels {^m rang- 
ing from 5 to 7 and a refinement criterion Aj ranging from 4 to 
10. As expected, at a constant Aj, the various AMR calculations 




2x105 3x10^ 
Porticles 



4x10* 



5x105 



Fig. 2. Ratio between the angular momentum 7(f) over the initial 
angular momentum Jo at different times for the SPH simulation 
with Ap = 5 X 10'' and An = 50. The ratio is plotted as a function 
of the number of particles, ordered in decreasing density. The 
value of the ratio is a mean value over 7 500 particles. The red 
curve coiTesponds to results at to. 



show a convergence. The SPH simulations were performed with 
a constant number of neighbors An - 50 and a total number of 
particles Ap ranging from 5 x 10^ to 5 x 10^. A consequence of the 
deteriorating resolution with increasing density in standard SPH, 
as mentioned earlier, is illustrated by the fact that, for equivalent 
initial condition, it is easier to get a better core resolution with 
AMR. 

In FiglHa) we show density profiles as a function of the ra- 
dius in the equatorial plane for AMR calculations with {^i„ = 7 
and Aj = 10 and for the SPH with Ap = 5x lO'^ and An = 50. The 
density profiles are very similar, indicating good convergence 
between the two methods. The behaviour differs at relatively 
high radius because of the external gas in the AMR method. In 
the present simulations, the dynamical time to reach is in- 
creased by ~ 5 kyr, because of the rotational support. As seen 
in Table |2] when one increases the resolution, one seems to con- 
verge toward this time. 

When the core becomes adiabatic, the angular momentum 
conservation induces the formation of an accretion disk around 
the central object. The centrifugal force becomes comparable 
to the gravitational one on the equatorial plane, slowing down 
the collapse. The outer collapsing gas, which has a supersonic 
infalling speed, meets suddenly the static gas of the core, cre- 
ating an accretion shock. This shock can be clearly seen in 
Fig- Eb) where the radial velocity component averaged over 
the equatorial plane is displayed. The accretion shock is de- 
scribed slightly more accurately with the AMR method. The 
SPH curve is smoother before the shock due to the artificial vis- 
cosity scheme. The slope before the shock strongly depends on 
the hydrodynamical solver so the results illustrate the difference 
between the hydrodynamical methods used in our two codes. 

4.2. Theoretical local angular momentum 

We now investigate the issue of angular momentum conserva- 
tion. Note that both SPH and AMR equations ensure conser- 
vation of the linear momentum. Considering our axisymmetric 
model, without azimuthal perturbation, we can easily investigate 
the effect of numerical resolution on angular momentum conser- 
vation. The local angular momentum should be well conserved, 
until azimuthal symmetry is broken. The loss of local angular 
momentum in our model is only due to unphysical transport in- 
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Fig. 3. Left plot: Cumulated hydrodynamical torque on the rotational axis for SPH calculations with A^p = 5 x 10^ and A^n = 50 
and at the same times as in Fig.|2] Right plot: Cumulated torque on the rotational axis due to the standard artificial viscosity for the 
same calculations and times. Particles are ordered in decreasing density and the torques are averaged over 7500 particles. At to the 
cumulated torques are significant for the densest particles. 



herent to the numerical methods used in the two codes. Thanks 
to its Lagrangian properties, the SPH calculations gives access 
to the angular momentum that each particle has initially, i.e. the 
angular momentum that particle should have if the numerical 
scheme was conserving it exactly. Having access to the parti- 
cle initial angular momentum, the loss of angular momentum 
is easily calculated. The azimuthal velocity component, directly 
linked to the angular momentum J, is given by 



Vfl 



xvy - yvx 



where r is the distance from the rotation axis 



and the angular momentum 



• yvx 



(12) 



(13) 



(14) 



The angular momentum conservation along the z-axis allows 
us to write 



(15) 



Thus, the theoretical angular velocity of a particle at time f 
is determined by the ratio between its initial angular momentum 
Jo = (xQVy o - yoVxfi) and its actual radius: 



Jo 

Vs.th - — ■ 
r 



(16) 



With Eq. ( fTSI l. we can compare the theoretical angular veloc- 
ity component to the numerical ones, and in particular with the 
AMR results for which we do not have access to a theoretical 
value. Since SPH and AMR density profiles are almost identical 
at to, we suppose that the previous mapping giving /o('"o) as a 
function of the radius in the SPH runs is also valid for the AMR 
calculations. Note also that the method cannot account for the 
displacement of the particles which would arise by changing the 
angular momentum. 



4.2.1. Azimuthal velocity component 

Figure|2]gives the ratio between the angular momentum J{t) and 
Jo for the SPH simulation with A^p = 5 x 10^ and An = 50. The 
particles are ordered in decreasing density and the ratio is aver- 
aged over 7500 particles. A first interesting result is that denser 
particles loose more angular momentum. At tff , denser particles 
have lost 3% of their initial angular momentum. In less resolved 
calculations, i.e. Np - 5x 10^ particles, the effect is stronger and 
the densest ones loose more than 10% of their initial angular mo- 
mentum in a free-fall time tff. This percentage slightly decreases 
when increasing the number of particles. The reason is that the 
denser the particle the larger the viscous torque (see below) and 
thus the larger the angular momentum transport. This numerical 
transport is amplified when the core is close to becoming adia- 
batic. 

Figure |3] shows the cumulated hydrodynamical (left-hand 
side) and viscous (right-hand side) torques on the rotational axis 
for the same SPH calculations at the same times as in Fig. |2] 
In principle, there should not be any torque on the rotational 
axis because of this axisymmetric model. The cumulated torques 
are computed and summed at each timestep for each particle. 
The value plotted is an average over 7 500 particles and par- 
ticles are ordered in decreasing density. Denser particles have 
the largest cumulated hydrodynamical and viscous torques. The 
friction forces corresponding to the hydrodynamical torque are 
stronger for these particles, which is due to the strong differential 
velocity. 

Figure |4^ displays the azimuthal velocity component as a 
function of the radius r on the xy-plane for the SPH simulations. 
The theoretical azimuthal velocity profile (solid line) is obtained 
following the previous section. It is obvious that low resolution 
simulations are not able to conserve properly the angular mo- 
mentum. With 5 X 10^ particles, we obtain counter-rotating par- 
ticles at the center (not illustrated in Fig. |4^ because of the av- 
erage in the xy-plane that smoothes the profiles). It appears that 
a minimum of 5 x lO'* particles is required to maintain angu- 
lar momentum loss within less than 10%, for the case of the 
present study. The improvement of angular momentum conser- 
vation eventually saturates for large numbers of particles. We 
checked that using a larger number of neighbors does not im- 
prove the conservation of local angular momentum. 
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Fig. 4. Azimuthal velocity at to as a function of the radius on the equatorial plane for SPH (left) and AMR (right) calculations 
at corresponding to. The left-hand plot (Fig. |4^) shows SPH results with various A^p and A^n - 50. The solid line represents the 
theoretical azimuthal velocity interpolated at to and is denoted as ve th- The right-hand plot (Fig. lit) shows AMR results with 
A^j - 10 and ^mi,, = 5, 6 and 7. The theoretical azimuthal velocity is plotted also for easy comparison with the SPH results. 



Table 2. Summary of the different simulations (upper table: 
SPH, lower table: AMR) performed to study angular momentum 
conservation. 





N, 


A'n 


Ni 


core 


to (kyr) 






5 X 10' 


50 


1.86 


225 


115 






1 X lO"* 


50 


2.34 


422 


107 






5 X lO'' 


50 


4. 


1 833 


98 






2x 10^ 


50 


6.35 


7 055 


95 






5x 10^ 


50 


8.61 


17 309 


93 


















/ - 

*^ mm 


A'i 


A^j 


^ core 


Tot. cells 


to (kyr) 


5 


2 145 


6 


3 928 


~ 9.1 X lO* 


150 


5 


2 145 


10 


30 752 


~ 1.6 


X 10^ 


116 


6 


17 160 


4 


4016 


~ 3.1 X 10' 


116 


6 


17 160 


10 


28 800 


~ 3.7 X 10^ 


109 


7 


137 260 


10 


29 944 


~ 2.2 X 10** 


96 



0.0010 r 



0.0001 



Np=5xl05, N„=50 

Lmn = 6, Nj=10 



14.0 14.5 15.0 15.5 16.0 16.5 

og(r) (cm) 



Figure ^ shows results obtained with the RAMSES code. 
AMR curves are plotted and compared with the theoretical one 
obtained previously for SPH. The simulations with {^i„ = 6 and 
^min = 7 are close to the theoretical curve. In both AMR simu- 
lations with A^j = 10, dense core resolution is higher than with 
SPH. For an initial resolution of ^min=5, the AMR scheme does 
produce some angular momentum lag. This can be due to the 
fact that with a poor initial resolution, the interpolation of the 
gravitational potential tends to convert gravitational energy into 
rotational energy. The outer gas does not alter the angular mo- 
mentum conservation for AMR calculations because of its tiny 
density. 

In Fig.|5^,b, we plot the integrated mass and the ratio of nu- 
merical over theoretical angular momentum for AMR calcula- 
tions with ^min - 6, A^j = 10 and for SPH calculations with 
A^p = 5 X 10^, An - 50. In Fig. [5^, most of the mass of the 
forming disk remains within a radius ~ 1 x 10'^ cm, i.e. where 
AMR and SPH have an opposite behaviour. In Fig. [Sj), angu- 
lar momentum is clearly transported to the outer regions with 
SPH, whereas it has a slight trend to be transported to the inner 
regions with the AMR. The overall angular momentum is well 
conserved in both calculations, but local properties seem to be 
affected by initial resolution. SPH resolution is much better at 



1.00 r 



Ni=5x105, Nn=50 - 

Uin = 6, Nj=10 



14.0 14.5 15.0 15.5 16.0 16.5 

og(r) (cm) 

Fig. 5. (a): Total integrated mass in the equatorial plane at to as a 
function of the radius for AMR calculations with /"min = 6, Aj = 
10 (dashed-Hne) and for SPH calculations with Ap = 5 x 10^, 
An = 50 (dotted line). (FiglSj^): Ratio between numerical and 
theoretical azimuthal velocities at to as a function of the radius 
in the equatorial plane for the same calculations. 



the beginning of the calculations, which enables SPH to prop- 
erly conserve local angular momentum at r > 1.4 x lO'*" cm at 
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Table 3. Summary of the different simulations performed with 
a - 0.5 and /? = 0.04. A^cme corresponds to particles/cells whose 
density satisfies p > 1 x 10"'^ g.cm"'' at t=to + 4 kyr. The upper 
table (Table |3^) gives a summary of the SPH calculations. The 
lower Table (Table[3j3) shows a summary of AMR calculations. 





A'n 


J ' core 


to (kyr) 


5x 10"* 


50 


6 997 


64 


2x 10^ 


50 


28 516 


63 


5x 10^ 


30 


72 417 


63 


5 X 10' 


50 


71 804 


62 


5 X 10' 


100 


76 390 


63 


5 X 10' 


200 


71 940 


64 



/ . 


A'j 


N 

^ ' core 


Tot. cells 


to (kyr) 


6 


4 


48 680 


~ 4.4 X 10' 


66 


6 


10 


156 588 


~ 6.2 x 10' 


65 


6 


15 


263 304 


~ 7.9 X 10' 


67 


7 


10 


117 108 


~ 2.4 X 10« 


65 


7 


15 


305 896 


~ 2.7 X 10* 


65 



to while low initial resolution of the sphere in AMR induces a 
worse conservation. This effect reverses at lower radius where 
AMR can reach smaller scales contrary to standard SPH. 

The density profiles obtained with the two methods converge 
towards a similar solution (c.f. FiglT^) but a closer analysis of 
the velocity profiles shows discrepancies. We can say that an- 
gular velocity profiles indicate that local angular momentum is 
better conserved with the AMR than with the SPH method for 
the code implementation that we used. Both can be improved 
using appropriate methods (see appendixiBlandlCli. 

5. Fragmentation 

5.1. Model 

Prestellar core fragmentation is a highly non-linear process. A 
key issue is to understand to what extent the fragmentation 
which occurs in numerical simulations is influenced by the nu- 
merical scheme and resolution. 

To study dense core fragmentation, we choose the same pre- 
vious spherical model and impose a m = 2 azimuthal density 
perturbation: 

p(6)) =po[l -HAcos(m0)], (17) 

where po is the mean sphere density, A the perturbation ampli- 
tude and 9 the azimuthal angle in cylindrical coordinates. 

The initial conditions are easy to implement for the AMR 
calculations. The SPH sinusoidal density perturbation is im- 
posed by adjusting the unperturbed 0-coordinate of each particle 
to a perturbed value 6* given by: 

A sin(m0*) 

0^0* + 1 i. (18) 

m 

We compare simulations for three different thermal supports, 
namely a = 0.35, 0.5 and 0.65, with a fixed rotational support 
j6 = 0.04. In the following, the convention is to call fragments 
the clumps where the gas density satisfies p > 1 x 10 g.cm 

5.2. Results for a critical case: a = 0.5,/? = 0.04 and A=0. 1 

This subsection is devoted to the exploration of various numer- 
ical parameters. First, we study the effect of varying the ini- 



tial grid resolution t^in and the number of cells within a Jeans 
length A^j for AMR calculations. Then, we present our SPH cal- 
culations with various number of neighbors A^n and of parti- 
cles A^p. For this set of calculations, the initial parameters are: 

po = 1.35x10-1'^ g.cm-^/^o = 7.07x10'^ cm, Oo = 2.12x10-'^ 
rad.s"' and tji = 57 kyr. The initial perturbation amplitude is 
A =0.1. 

Table [3] summarizes the calculations we performed for this 
case. Informations about the core resolution (i.e. p > 1 x 10"'^ 
g.cm"^) and the total number of cells are given at t=to+4 kyr. The 
dynamical times to reach the collapse are quite similar, within 
less than 2%. Synchronizing calculations at to is then well justi- 
fied. In the following sections, we consider core evolutions over 
a few thousand years (~ 10% of to). 

5.2.1 . Detailed study of the effect of A^j and 4,in for AMR 
calculations 

In Fig.|6] we show density maps on the equatorial plane at four 
different timesteps for AMR calculations with, from left to right, 
^min = 6 and Nj = 4, 10 and 15 and Anin = 7 (giving 128"* cells 
initially) and A^j = 10. Maps are given, from top to bottom, at 
t=to + 4 kyr, t=to + 5 kyr, t=to + 6 kyr, and t=to + 7 kyr. As shown 
in Fig. |6] fulfilling the Truelove condition (Nj > 4) does not 
guarantee an accurate fragmentation timescale. The calculations 
with Aj = 4 fragments at to + 5 kyr while other calculations do 
not fragment until t=to-i-7 kyr This suggests that calculation with 
fmin = 6 and Aj = 4 suffers from inaccurate fragmentation, but 
this latter is inhibited when fragments fall on the central object 
before to + 6 kyr The core will eventually refragment but not 
at the same time as the other calculations (i.e. t >> to -i- 7 kyr). 
With increasing Nj, we converge to a fragmented pattern with 
one central object and two satellites. 

Another aspect to be considered quite carefully is the choice 
of imin, i c- the initial description of the sphere. The two calcu- 
lations with Aj = 15 and ^min = 6 and £„iin - 7 are very similar, 
suggesting that numerical convergence has been achieved. Even 
though small differences still appear in the detailed structures. 
The satellites formed with fn,in = 7 are more structured and 
compact than those formed with the initial resolution t^nm = 6. 
According to these calculations, fragmentation into two identi- 
cal satellites and a central object should occur around t=to + 7 
kyr. 

5.2.2. Detailed study of the effect of Ap and An for SPH 
calculations 

We performed a series of SPH calculations with different Ap, but 
with a constant An = 50, and then with a larger An for one value 
ofAp. 

Figure Q shows density maps on the equatorial plane for 
these calculations at four different timesteps, namely, from top 
to bottom: t=to + 4 kyr, t=to + 5 kyr, t=to + 6 kyr and t=to + 7 kyr 
The calculations with Ap = 2x10^ and Ap =5x10^ show good 
agreement at least until to + 6 kyr and differ noticeably from the 
less resolved calculations (An = 50, Np - 5 x 10'*) where frag- 
menatation occurs earlier. The dense core fragments in a ny case, 
but its fragmentation is delayed when Ap increases (e.g. NelsonI 
|2006). Early fragmentation is here clearly due to a lack of resolu- 
tion. As shown in section |4.2.1| the conservation of local angular 
momentum is bad when Ap is low and leads to very inaccurate 
collapse and fragmentation timescales of the cloud. Once sym- 
metry is broken, it is useless to continue the simulations, since 
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Fig. 6. AMR calculations density maps on the xy-plane for the case a - 0.50, /3 = 0.04. From top to bottom, four different times are 
showed: t=to + 4 kyr, t=to + 5 kyr, t=to + 6 kyr and t=to + 7 kyr. The AMR calculations have been performed with, from left to right 
columns, /"min - 6 and A^j = 4, A^j = 10, Nj = 15 and ^min - 7 and Nj - 15. 



calculations would obviously diverge. This symmetry breaking 
occurs earlier in the SPH calculations first because of the nu- 
merical noise inherent to the relaxed, and random, initial parti- 
cle distributions, and also because of the lower resolution of the 
Jeans length in the disk. 



The other fundamental parameter in SPH calculations is the 
number of neighbors determining the kernel size. Increasing the 
number of SPH particles increases the resolution but also in- 
troduce numerical noise at smaller scales. The natural way to 
reduce this noise is to increase the smoothing kernel length by 
in creasing An. The effect s of varying An have been invest igated 
bvlLombardi et all d 19991) and lRasiol (ll999l) . In particular. [Rasiol 
derived the following results: 



- higher accuracy is reached when both Ap and An are 
increased, with Ap increasing faster than An so that 
the smoothing leng th decreases. One possible scaling 
dLombardi et al.|[T99 9) is An « A* with 0.2 < ^ < 1, 

- SPH scheme is consistent in the limit where (An,A^p) —> °° 
and h 0, 

- convergence (e.g. the number of timesteps) is accelerated by 
increasing the smoothness of the kernel. 

The usual number of neighbors in previous studies is about 
50. However, no study has really explored the role played by 
An in the context of star formation. Hence, we performed cal- 
culations with a constant Ap and different values of An. This is 
illustrated in the right column in Fig. |7] where we report maps 
of calculations with Ap = 5 x 10^ and An = 100. The core has 
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Fig. 7. SPH calculations density maps in the jcy-plane for the case a - 0.50, p - 0.04. From top to bottom, four different times are 
showed: t=to + 4 kyr, t=to + 5 kyr, t=to + 6 kyr and t=t() + 7 kyr The calculations have been performed with, from left to right, 
A^N = 50 and A^p = 5 x 10^ Ap = 2 x 10^ and Ap = 5 x 10^ and An = 100 and Ap = 5 x 10^ 



already fragmented at to + 6 kyr whereas with An — 50 it frag- 
ments later. The calculations with parameter set Ap = 2 x 10^ 
and An = 50 is very similar to the later of similar ratio Ap/AN. 
Other calculations with various An are reported in Appendix lAl 
It appears clearly that the greater An the earlier fragmentation 
occurs, because increasing An for a fixed Ap decreases the spa- 
tial resolution {h increases). 



5.2.3. Comparison and convergence 

In the previous sections, we show that AMR and SPH calcu- 
lations converged separately. We now cross-compare the con- 
verged calculations. Figure[8]shows density maps in the equato- 
rial plane for the results of two amongst the most resolved calcu- 



lations at three timesteps, namely, from top to bottom, to +5 kyr, 
to+6 kyr and to+7 kyr. The left column shows maps for AMR 
calculations with fn,in - 7 and Aj = 15 whereas the right col- 
umn displays SPH maps for calculations with Ap = 5 x 10^ and 
An = 50. We display again the results of Fig|6]and Fig|2]to illus- 
trate clearly the convergence. Agreement between the two meth- 
ods for these physical and numerical parameters set is striking 
for the two first timesteps. The calculations give the same frag- 
mentation time and pattern, although satellites and the central 
object are bigger with the SPH. 

Figure |9] shows disk density profiles as a function of the 
radius averaged in the equatorial plane for the same SPH and 
AMR calculations. This plot complements Fig. [8] with the last 
density maps, where fragments are well developed. Density 
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profiles show a peak at a radius corresponding to satellite 
positions in the map. Satellite fragments are denser in the SPH 
calculations, and the central object is less dense and bigger 
compared with the dense elongated shape obtained with the 
AMR. 

Although there are some obvious differences between the 
two methods , there seems to be a real convergence between the 
two types of calculations. For the specific case under study, we 
find a good agreement between AMR calculations with = 7 
or 6 and A^j = 15 and SPH calculations with A^p = 5 x 10^ and 
A^N = 50. However, even for the most resolved simulations, the 
results between the two methods diverge after some time (i.e. 
t{)+7 kyr for this specific case). This is not very surprising be- 
cause the dynamics becomes very non-linear and chaotic, the 
initial and numerical noise are getting amplified. 

Ni=5xl05, Nn=50 
t=67.01 kyr 



1883 _ Nj=15 
t=69.68 kyr 



P (g-cm 



t=70.77 kyr 



t=6B.17 kyr 
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X (AU) X (AU) 

Fig. 8. AMR and SPH calculations density maps in the xy-plane 
at three different times for the case a = 0.50, j6 - 0.04. The times 
correspond to to -n 5 kyr, to + 6 kyr and to -n 7 kyr, from to bottom, 
respectively. The AMR calculations plotted on the left column 
have been run with ^^in = 7 and Nj = 15. The right column 
shows the results for the SPH calculations with A^p = 5 x 10^ and 
An = 50. 



5.3. Results for low and high thermal support 

5.3.1 . Results for a least prone to fragment case: a = 0.65, 
yS = 0.04 

This second series of calculations is the least prone to fragmenta- 
tion because of its strong thermal support. We use a perturbation 
amplitude A - 0.5 in order to make fragmentation easier if it 
should occur 
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Fig. 9. Density profiles at to+ 7 kyr as a function of the radius, 
averaged on the equatorial plane, for AMR calculations with 
^min = 7 and Aj = 15 (dashed line), and SPH calculations with 
Ap = 5 X 10-** and An = 50 (dotted fine). 



Figure [TO] gives density slices on the equatorial plane at 
t=to -H 10 kyr On the left-hand side, we show AMR results for 
an initial sphere described with ^„„„ - 6 and, from top to bot- 
tom, Aj -4-, 10 and 15. The right-hand column shows slices of 
SPH calculations with An = 50 and Ap = 5 x 10"*, 2 x 10^ and 
5 X 10^, from top to bottom. In the previous cases, the core has 
akeady fragmented into three clumps at this time. In this case, 
the cloud develops spiral arms with no fragmentation. The cloud 
fragments in some cases after to + 50 kyr. 

AMR and SPH calculations converge quickly to a pattern 
with only spiral arms and the formation of a single central object 
when resolution is sufficient. 



5.3.2. Early fragmentation case: a - 0.35 , /? - 0.04 

This last case is the most prone to fragmentation because of its 
small thermal support against gravitational energy. In this set 
of calculations, the initial parameters are: po = 3.92 x 10"'^ 
g.cm-^ Ro = 4.95 x 10"' cm, Oq = 3.63 x 10"'^ rad.s"' and 
tff = 1.06 X 10'^ s (~ 33.6 kyr). The initial perturbation ampli- 
tude is A =0.1. We performed SPH calculations with Ap ranging 
from 5 X 10"* to 5 X 10^ and An = 50. The AMR calculations were 
performed with ^mh, - 6 and 7 and Aj varying between 4 and 15. 
Although all results are very similar at to, we find some differ- 
ences at to -I- 1 kyr For example, it appears quite clearly that 
AMR calculations with - 6 and Aj = 4 diverge from the 
other AMR calculations (tiny spiral arms). 

Figure [TT] shows density maps in the equatorial plane for the 
most relevant calculations at t=to + 3 kyr On the left column, 
we give the AMR results with increasing resolution parameter 
Aj from top to bottom and a constant {^in = 6. According to our 
previous results, an initial computational domain with {^in = 6 
is sufficient to reach convergence for the AMR calculations. The 
right column shows SPH calculations, with Ap ranging from 
5 X 10"* to 5 X 10^. We seem to reach a convergence between 



the AMR calculations with 



= 6 and Aj > 6. The SPH cal- 



culations with Ap = 5 X lO'' diverge quickly compared to two 
more resolved with Ap - xlO^ and 5 x 10^'. The most resolved 
AMR and SPH runs show a convergence towards a similar so- 
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Fig. 10. Density maps in the equatorial plane at to + 10 kyr for Fig. 11. Density maps in the equatorial plane at t=to + 3 kyr for 
a - 0.65,13 - 0.04 and A = 0.5. On the left-hand side, we show the case a - 0.35, jS - 0.04 and A - 0.1. On the left side. 



AMR results with 4 



6 and, from top to bottom, N] - 4, we give AMR results with 



6 and, from top to bottom. 



10 and 15. The right-hand side gives slices of SPH calculations N] = 4, 10 and 15. On the right side, we show SPH maps with 
with A^N - 50 and A^p = 5 x 10"*, 2 x lO*" and 5x10^, from top calculations with An = 50 and, from top to bottom, Ap =5x10'*, 



to bottom. 



An = 2 X 10-** and An = 5 x 10^. 



lution. The patterns have the same size and position. The core 
fragments into a central clump (of size ~ 30 AU) and two iden- 
tical outlying clumps (of size ~ 10 AU) for AMR calculations. 
SPH results give a similar central object, but the outlying clumps 
are larger 

However, as shown in Appendix C, higher resolution runs 
show that convergence has not been reached. In that case, one 
needs either an even better resolution or, alternatively, a more 
powerful numerical scheme. 

To conclude our study on core fragmentation, we can say 
that the more non linear is the issue, the more difficult it is to get 
convergence between SPH and AMR simulations. Good conver- 
gence is found for high enough thermal support. However, for 
low thermal support (i.e. a - 0.35), convergence is more diffi- 
cult to achieve. The horizon of predictability in such a case is 
very short. 

6. Summary and Discussion 

We have investigated the effect of numerical resolution in AMR 
and SPH calculations on the collapse and the fragmentation of 
rotating cores. 



We show that we reach good convergence between AMR and 
SPH methods provided one uses sufficient numerical resources. 
First, we take a simple model to study local angular momen- 
tum conservation. The initial study shows that local angular mo- 
mentum is better conserved with the AMR approach for equiv- 
alent computational needs, whereas SPH gives better dynami- 
cal times. As shown in Fig. |4^, a smaller number of particles 
in standard SPH calculations leads to bad local angular momen- 
tum conservation. Numerical torques on the rotational axis are 
accumulated for denser particles, whereas our model should re- 
main axisymmetric. In AMR calculations, a poor initial com- 
putational domain resolution (i.e. i^i^ < 6) leads to unphysical 
transfer of gravitational energy to rotational energy (see Fig.|4j3). 
A significant loss of angular momentum will affect fragmenta- 
tion since less rotational support can balance gravitational col- 
lapse. The smallest parameter set for SPH calculations required 
to go through gravitational collapse without significant loss of 
angular momentum corresponds to a number of ~ 530 particles 
per Jeans mass at the critical density pc, i-e. 5 particles per Jeans 
length. The equivalent minimum resolution criterion for AMR 
calculations is ^,ni„ 6 and Aj - 4. 
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Then we investigate fragmentation issues for three different 
initial condition. For the least prone to fragment case (a - 0.65, 
y6 - 0.04), we show that AMR and SPH methods give similar 
results when the respective Jeans resolution criteria are fulfilled. 
These results agr ee with the semi-analytical c riteria on fragmen- 
tation derived by iTsuribe & Inutsukal (1 19991) for the isothermal 
collapse (a ^ 0.55 - 0.65 yS). We study extensively the case 
a = 0.5,/? = 0.04. We first reach good agreement between AMR 
calculations with the parameters ^min - 6 and A^j = 15 and SPH 
calculations with A^p = 5 x 10^ and An = 50, i.e. ~ 5370 par- 
ticles per Jeans mass at critical density pc. These parameter sets 
seem to be a lower resolution limit for dense core collapse and 
fragmentation SPH and AMR calculations in order to get good 
agreement in both time and space scales. Using a lower number 
of particles or number of points per Jeans length will lead to in- 
accurate early fragmentation due to numerical effects. Initially, 
we compare the two converged calculations and, for this specific 
case, we find good agreement between the two methods (see Fig 
[8]). The price to pay in computer time, however, is larger with 
the SPH method for the fragmentation study, due to not using 
sink particles. In the case of low thermal support, the dynamic 
quickly becomes very non-linear and numerical convergence of 
the simulation can not be achieved as easily as for higher thermal 
support. A statistical analysis over a large number of simulations 
would be needed to see if converging results can be obtained in 
term of fragment distributions (in mass, size ...). 

The two approaches show good agreement for the gen- 
eral pictures. Details are better resolved in AMR calculations 
thanks to the refinement method based on the local Jeans length, 
whereas the resolution deteriorates with increasing density with 
standard SPH. Numerical calculations of protostellar collapse 
should thus be conducted with great care, with a detailed exami- 
nation of numerical resolution. The present work can be used to 
assess the validity of numerical tools to study star formation. 
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Appendix A: Complementary results on the effect 
of for SPH calculations for the case: 

a - 0.5,/3 - 0.04 and A=0.1 

Figure lA.ll shows SPH calculations run with a constant Np = 
5x10^ and values of An = 30, 50, 100 and 200 from top to bot- 
tom for two timesteps (t()-H5 kyr on the left column and t()-i-6 
kyr on the right column). These simulations should be com- 
pared with simulations presented in Fig. |2l The first relevant 
result is the fact that increasing An speeds up fragmentation. 
Moreover, there seems to be a similarity between calculations 
with low Ap/AN ratio, i.e. An = 50, Ap = 2 x 10^ on one side 
and An = 100 and Ap = 5 x 10^ on the other side. We find 
the same patterns at different times, postponed when either Ap 
increases or An decreases, the number of resolution elements 
being equal. This illustrates the compromise between resolu- 
tion and convergence that must be respected in SPH calculations 
(Lombardi et al. . 1999: .Rasio. 1999,) . 

Appendix B: Note on the artificial viscosity and 
numerical diffusion in SPH 

Diffusivity is a well-known drawback of standard SPH. This is- 
sue can be reduced us ing a constant number neighbors, AAn - 
dAttwood et al.ll2007l) . and ad vanced scheme for viscosity such 
as time-dependent viscosity dMorris & MonaghanI 11997 ). We 
present here SPH calculations of the collapse of the uniform- 
density sphere already studied in ^ but using another scheme 
viscosity and/or a constant number of neighbors. This two im- 
provements are quite easy to implement and do not require 
expensive extra computational costs. Time-dependent viscosity 
calculations have been done with a* = 0.1, and an e-folding 
constant equal to 0.15. 

Figures IB. lb shows the averaged ratio between angular mo- 
mentum J{t) at time to and initial angular momentum Jo as a 
function of particles (ordered in decreasing density) for SPH 
calculations run with Ap =5x10^ and An - 50 and the im- 
provement above mentioned either turned on or not. It is clear 
that time-dependent viscosity better conserve angular momen- 
tum for the denser particles. Hence, less angular momentum is 
transported to the outer part of the core. This is confirmed in Fig. 
IB. lb where we plot the ratio between angular momentum at time 
to and initial angular momentum as a function of the radius for 
the same calculations. Local angular momentum conservation is 
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Fig.A.l. Density maps in the equatorial plane at two different 
times from SPH calculations with a - 0.50, yS = 0.04 and 
A^p = 5 X 10^. The left column shows density maps for calcu- 
lations with A^N - 30, 50, 100 and 200 from top to bottom, re- 
spectively, at t=to + 5 kyr. The right column represents the same 
calculations at t=to -i- 6 kyr. 



increased by 10% in the inner part. However, keeping constant 
the number of neighbors does not improve local angular momen- 
tum conservation sinc e the system only evo lves over about one 
free-fall time whereas [Attwood et al.l (|2007|) shows that dissipa- 
tion becomes significant after a few free-fall times. All these im- 
provements of standard SPH are as many new features that will 
strengthen convergence with the AMR, particle-spUtting being 
the most promising one. 
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Fig.B.l. Top plot (Fig IrH i): same as Fig. |2] with different AAn 
and artificial viscosity scheme, yiif) indicates the use of time- 
dependent viscosity instead of the standard artificial viscosity 
scheme in the paper Bottom plot (Fig IB. lb ): same as Fig.|5]for 
the above mentioned SPH calculations. 



Appendix C: Note on the diffusion of the numerical 
schemes in AMR 

A key ingredient in the AMR method is the numerical scheme 
used to compute flux at the grid's interfaces. In this paper, we 
use a Lax-Friedrich ( hereafter LF) Rieman n solver designed 
for MHD calculations (iFromang et al.ll2006l) . However, the LF 
scheme is known to be a diffusive scheme. In this appendix, we 
present AMR calculations carried out with, on one hand, the LF 
scheme and, on the other hand, a Roe scheme. Roe scheme be- 
ing less diffusive than LF, this could have dramatic effect on the 
fragmentation issue. 



C.7. Case a = 0.5, /? = 0.04 

Figure ICTI shows density maps on the equatorial at three differ- 
ent timesteps for two AMR calculations run with the same nu- 
merical parameters, i.e. ^„„„ = 6andAj = 15, but with a different 
solver, i.e. the LF one on the left column and the Roe one on the 
right column. Results are quite similar, AMR calculations are in 
good agreement for this critical case with the two solvers. Since 
less angular momentum has been locally lost or transported with 
the Roe scheme, the core is smaller and the fragments are closer 
to the central object. This brings support to the fact that we find 
good convergence between AMR and SPH calculations for this 
case. 
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Fig. C.l. Density maps in the equatorial plane for the case a - 
0.5, p - 0.04 at, from top to bottom, t=t() + 5 kyr, t=to + 6 kyr 
and t=to + 7 kyr On the left-hand side, AMR results with Lax- 
Friedrich solver and tmin = 6, A^j = 15 are reported and AMR 
results of calculations with he same parameters but with a Roe 
solver are given on the right column. 



C.2. Case a = 0.35, yS = 0.04 



Figure IC.2I shows density maps on the equatorial plane at to + 
2 kyr (right column) and to + 3 kyr for three simulations of the 
case a - 0.35, p - 0.04 with numerical parameters ^min - 6 
and Nj = 12 (top and bottom maps, LF and Roe schemes) and 
^min = 7 and A^j = 15 (middle row, LF scheme). Let us remind 
that in Fig. (TT] the case £mm - 6 and Aj = 15 with the LF 
solver has been displayed. The two calculations with the same 
numerical parameters differ, according to the numerical scheme 
used. The fragmentation process changes: one gets a configura- 
tion central object + two satellites with the LF scheme whereas 
we get a binary system resulting from the fragmentation of the 
central object with the Roe scheme. If we improve the initial 
sphere resolution in LF calculations (i.e. {^m = 7, Aj = 15), 
we converge to the results obtained with the Roe scheme, i.e. 
a central binary system, with a value ^min < 7. We know that 
angular momentum is well conserved using the Roe scheme or 
the LF scheme with ^,nin - 7, so it seems that calculations lead 
to a different core fragmentation because of their less accurate 
angular momentum conservation. Since we use a small thermal 
support, it is easy to reach another fragmentation configuration, 
these processes being highly non-linear. 



Appendix D: Note on SPH sink particles 

The introduction of sink particles is a widely used way to get 
a compromise between good resolution and acceptable timestep 



IB-'* LD-" JD-"* m-"' ID-"* ID--" in-"J 10-"' lo-" 



♦ 



<f\<b 



R., \f ja, (tyr, I n. 0:^73 



B, Vr [a. Unr:, I fl.flWJ 



-qw -IK 4 l« 



-lit t Liyt ^> 



Fig. C.2. Density maps in the equatorial plane for the case a = 
0.35, p - 0.04 at t=to + 2 kyr on the left-hand side and t=to + 
3 kyr on the right-hand side. For the two upper rows, we plot 
results for AMR calculations with (£mm = 6, Aj = 12), (^^in = 
7, Aj = 15) and our usual Lax-Friedrich scheme. The bottom 
row gives the results for calculations conducted with (fmin = 6, 
Aj = 12) too, but with a Roe solver. Times are given in Myr. 



in SPH methods. Creating a sink particle enables to loosen the 
Courant condition on the particle timesteps. 

The density level psink at which a sink particle is created has 
to be chosen with care. In the previous SPH calculations, no sink 
particles were used. Let us focus on the highly non-linear frag- 
mentation case, a = 0.35, y6 = 0.04 and A = 0.1 to present SPH 
calculations carried out with various sink densities. 

Figure lDTI shows calculations carried out with three different 
densities for the creation of sink particles, psink, namely 1 x 10"'° 
g.cm"-' (resulting in no sink creation), 1 x 10"" g.cm"-' and 3 x 
10"'^ g.cm"'-, from top to bottom, and the same number of SPH 
particles and neighbors, Ap = 5 x 10^ and An = 50. The left- 
hand side reports results at to-H 1 kyr and the right-hand column 
results at to H- 3 kyr The CPU time is about 28% smaller for the 
case Psink = 3 X 10"'^ g.cm"-' than for the case psink = 1 x 10""^ 
g.cm"^, but the dense core resulting at to+3 kyr is really different. 
The upper row corresponds to calculations without creation of 
sink particles, whereas one sink particle has been created with 
Psink - 1 X 10"" g.cm"^ and 9 with psi„k = 3 x 10"'^ g.cm""*. It 
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is easy to see that, even if only one sink particle is created, the 
complete dynamic is affected particularly in the central region. 

p (g.cm"^) p (g.cm"^) 
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Fig. D.l. Density maps in the equatorial plane for the case a = 
0.35, /? = 0.04 at t^to + 1 kyr on the left-hand side and t=to -H 
3 kyr on the right-hand side. 9 sink particles have been created 
on the bottom figures, affecting the whole dynamic of the dense 
core. 

In conclusion, it is clear that sink particles should be handled 
with great care. A fair comparison should also compare these 
SPH calculations with AMR ones including sink particles. Such 
studies are in progress. 



